The objective of this example is to demonstrate how to use the catalogue tools to calculated recurrence by different methodologies
In [1]:
# Import required modules
import numpy as np
import matplotlib.pyplot as plt
from hmtk.parsers.catalogue.csv_catalogue_parser import CsvCatalogueParser, CsvCatalogueWriter
# Import Recurrence Tools
from hmtk.seismicity.occurrence.b_maximum_likelihood import BMaxLikelihood
from hmtk.seismicity.occurrence.kijko_smit import KijkoSmit
from hmtk.seismicity.occurrence.weichert import Weichert
# Import Mmax Tools
from hmtk.seismicity.max_magnitude.kijko_nonparametric_gaussian import KijkoNonParametricGaussian
from hmtk.seismicity.max_magnitude.kijko_sellevol_bayes import KijkoSellevolBayes
from hmtk.seismicity.max_magnitude.kijko_sellevol_fixed_b import KijkoSellevolFixedb
from hmtk.seismicity.max_magnitude.cumulative_moment_release import CumulativeMoment
print 'Import OK!'
In the present example we use a synthetic catalogue with known b-value properties and completeness windows
In [3]:
input_file = '../../hmtk/tests/seismicity/completeness/data/completeness_test_cat.csv'
parser = CsvCatalogueParser(input_file)
catalogue = parser.read_file()
print 'Input complete: %s events in catalogue' % catalogue.get_number_events()
print 'Catalogue Covers the Period: %s to %s' % (catalogue.start_year, catalogue.end_year)
# Plot magnitude time density
from hmtk.plotting.seismicity.catalogue_plots import plot_magnitude_time_density
magnitude_bin = 0.1
time_bin = 1.0
plot_magnitude_time_density(catalogue, 0.1, 1.0)
In [4]:
# Known completeness properties
completeness = np.array([[1990., 4.0],
[1960., 4.5],
[1910., 5.5]])
from hmtk.plotting.seismicity.catalogue_plots import (plot_observed_recurrence,
get_completeness_adjusted_table,
_get_catalogue_bin_limits)
earthquake_count = get_completeness_adjusted_table(catalogue, completeness, 0.1, 2009)
print 'Magnitude N(OBS) N(CUM) Log10(Nc)'
for row in earthquake_count:
print '%6.2f %10.3f %10.3f %10.3f' %(row[0], row[1], row[2], row[3])
plot_observed_recurrence(catalogue, completeness, 0.1)
The first method is simply an adaptation of the Aki (1965) approach in which b-value is determined from the weighted mean of the b-value for each completeness period.
In [5]:
# Set up the configuration parameters
recurrence_config = {'reference_magnitude': None,
'magnitude_interval': 0.1,
'Average Type': 'Weighted'}
bml_recurrence = BMaxLikelihood()
bval, sigmab, rate, sigma_rate = bml_recurrence.calculate(catalogue, recurrence_config, completeness)
print 'B-value = %9.4f +/- %9.4f' %(bval, sigmab)
print 'Rate (M >= 4.0) = %9.4f +/- %9.4f' %(rate, sigma_rate)
Implements the maximum likelihood estimator of Kijko & Smit (2012)
In [6]:
# Set up the configuration parameters
bks_recurrence = KijkoSmit()
bval, sigmab, rate, sigma_rate = bks_recurrence.calculate(catalogue, recurrence_config, completeness)
print 'B-value = %9.4f +/- %9.4f' % (bval, sigmab)
print 'Rate (M >= 4.0) = %9.4f +/- %9.4f' % (rate, sigma_rate)
Implements the maximum likelihood estimator of Weichert (1980)
In [7]:
# Set up the configuration parameters
bwc_recurrence = Weichert()
bval, sigmab, rate, sigma_rate = bwc_recurrence.calculate(catalogue, recurrence_config, completeness)
print 'B-value = %9.4f +/- %9.4f' %(bval, sigmab)
print 'Rate (M >= 4.0) = %9.4f +/- %9.4f' %(rate, sigma_rate)
In [ ]:
Maximum likelihood estimator of maximum magnitude assuming no uncertainty in b-value (see Kijko (2004) for more details)
In [8]:
mmax_config = {'b-value': 1.0,
'input_mmin': 4.5,
'input_mmax': None,
'input_mmax_uncertainty': None}
mmax_ks = KijkoSellevolFixedb()
mmax, mmax_sigma = mmax_ks.get_mmax(catalogue, mmax_config)
print 'Mmax = %8.3f +/- %8.3f' %(mmax, mmax_sigma)
Maximum likelihood estimator of maximum magnitude assuming uncertain b-value (see Kijko (2004) for more details)
In [9]:
mmax_config = {'b-value': 1.0,
'sigma-b': 0.05,
'input_mmin': 4.5,
'input_mmax': None,
'input_mmax_uncertainty': None}
mmax_ksb = KijkoSellevolBayes()
mmax, mmax_sigma = mmax_ksb.get_mmax(catalogue, mmax_config)
print 'Mmax = %8.3f +/- %8.3f' %(mmax, mmax_sigma)
Maximum likelihood estimator of maximum magnitude assuming no specified magnitude frequency distribution (see Kijko (2004) for more details)
In [10]:
mmax_config = {'number_earthquakes': 100, # Selects the N largest earthquakes in the catalogue for analysis
'input_mmax': None,
'input_mmax_uncertainty': None}
mmax_knpg = KijkoNonParametricGaussian()
mmax, mmax_sigma = mmax_knpg.get_mmax(catalogue, mmax_config)
print 'Mmax = %8.3f +/- %8.3f' %(mmax, mmax_sigma)
Estimator of Maximum Magnitude based on the "Cumulative Moment" method, an adaptation of the cumulative strain energy estimator of Mmax (Makropoulos & Burton, 1983), with Mmax uncertainty estimated via Monte Carlo sampling of the observed magnitude uncertainties
In [11]:
mmax_config = {'number_bootstraps': 1000} # Number of samples for the uncertainty analyis
mmax_cum_mo = CumulativeMoment()
mmax, mmax_sigma = mmax_cum_mo.get_mmax(catalogue, mmax_config)
print 'Mmax = %8.3f +/- %8.3f' %(mmax, mmax_sigma)